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Abstract 

We investigate the Equation of State (EOS) of classical systems having 300 
and 512 particles confined in a box with periodic boundary conditions. We 
show that such a system, independently on the number of particles inves- 
tigated, has a critical density of about 1/3 the ground state density and a 
critical temperature of about 2.5 MeV . The mass distribution at the critical 
point exhibits a power law with r = 2.23. Making use of the grand partition 
function of Fisher's droplet model, we obtain an analytical EOS around the 
critical point in good agreement with the one extracted from the numerical 
simulations. 



PACS: 25.70.-z, 25.75.+r, 64.30. +t, 64.70.-p 



*present address: Universitat Tubingen, Institut fur Theoretische Physik, D-71076 Tubingen, 
Germany 



1 



I. INTRODUCTION 



Recently, several experiments in proton-nucleus and nucleus-nucleus reactions have re- 
vealed the creation of fragment size distributions exhibiting power laws, expected according 
to the droplet model of Fisher, for fragment formation near the critical point of a liquid- 
gas phase transition M. This observation raised the question whether critical phenomena, 
which strictly speaking can only occur in the thermodynamic limit, can be explored within 
the context of heavy- ion collisions, where only a few hundreds of nucleons are involved M . 

In previous works by studying the expansion of a classical finite system within the 
framework of Classical Molecular Dynamics (CMD), evidence has been found for the oc- 
currence of a critical behaviour of our system revealed through a study of fragment size 
distributions, scaled factorial moments and moments of size distributions. Such a critical 
behaviour is connected to a liquid-gas phase transition by the use of Fisher's droplet model 
Q and Campi analysis ||. In the present work, using the same two-body interaction ||, 
we calculate by means of the virial theorem |2||| , the equation of state of a classical system 
confined in a cubic box with constant density and periodic boundary conditions. This equa- 
tion of state has a critical temperature of about 2.5 MeV and a critical density of about one 
third the ground state density. Trajectories in the (p, T) plane of expanding finite systems 
indeed go through this critical point 0. It is also shown that at the critical point, the 
fragment size distribution exhibits a power law with r = 2.23, consistent with the droplet 
model of Fisher ||,|J, and scaling laws of critical exponents for a liquid-gas phase transition. 

We calculate in Section II, within the framework of CMD |||6|], the equation of state of 
the system in a cubic box with cyclic boundary conditions. Section III contains the study of 
fragment size distributions at the critical point and around it. In Section IV, we study the 
critical behaviour of the moments of size distributions using the grand partition function of 
the Fisher's droplet model and extract the equation of state of this model in the vicinity of 
the critical point, which is found in good agreement with the one calculated for the infinite 
system. Finally, we give our summary and conclusions in Section V. 
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II. EQUATION OF STATE 



In the following, making use of the virial theorem P||, we study the equation of state 
of classical systems with 300 particles (150 protons + 150 neutrons) and 512 particles (256 
protons + 256 neutrons) confined in a cubic box with periodic boundary conditions. The 
particles of the system move under the influence of a two-body potential V given by || : 

V np (r) = V r [exp(-fi r r)/r - exp(-fi r r c )/r c ] 

- V a [exp(-fi a r)/r - exp(-fi a r a )/r a ] 
V nn (r) = V pp (r) = V [exp(-fi r)/r - exp(-fi r c )/r c ] (1) 

r c = 5.4 fm is a cutoff radius. V np is the potential acting between a neutron and a proton 
while V nn is the potential acting between two identical nucleons. The first potential is 
attractive at large r and repulsive at small r, while the latter is purely repulsive so no bound 
state of identical nucleons can exist. The various parameters entering eq.(l) are defined 
with their respective values in Ref. ||. The classical Hamilton's equations of motion are 
solved using the Taylor method at the order 0[(5t) 2 ] where 5t is the integration time step 
M. Energy and momentum are well conserved because we use a very small time step St 
varying in time according to the maximum acceleration of the particles. 

Two different methods were used for initialization and thermalization of the system. In 
the first one, the particles are initially distributed randomly in a cubic box with a small 
size (the initial density is about 0.32 fm~ 3 ) and they are given the same velocity modulus 
which is quite large. The system is then cooled down to a predefined average kinetic energy 
by rescaling the velocities at each time iteration by a coefficient close to unity. After the 
density is decreased slowly to the desired density (at each time iteration the size of the box 
is slowly increased). During the cooling process, the system also thermalizes. Nevertheless, 
after the desired density and temperature have been reached, we let it evolve for additional 
time to allow a complete thermalization. In the second method, the particles are initially 
distributed on a cubic lattice such as the nearest neighbors of a proton are all neutrons 
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and vice versa, and the system is excited at an initial temperature giving to its particles 



then evolves according to Hamilton's equations of motion and reaches equilibrium after few 
thousands of iterations. With the first method, the system reaches equilibrium faster than 
with the second one. 

In Fig. 1 we have plotted the time evolution of the kinetic energy (upper part), the 
potential energy (lower part, solid line) and total energy (lower part, dashed line) computed 
using the second method for thermalization for the system with 512 particles. The density of 
the system is fixed at p = 0.16 fm~ 3 and the initial temperature is set at T = 1.66 MeV. We 
see that the total energy is well conserved during all the time evolution (2.5% in 100 fm/c). 
In this case the system thermalizes very quickly in the first 20 fm/c, then the kinetic energy 
stabilizes around 1.0 MeV (T = 0.7 MeV) with small fluctuations due to the interplay 
between kinetic and potential energy. In Fig. 2 we plot the same quantities at a density of 
0.04 fm~ 3 and the same initial temperature T = 1.66 MeV. In this case we see that the 
thermalization process takes a longer time and the kinetic energy stabilizes around 4. MeV. 

In the following, the calculations are done using the first method of initialization and 
thermalization for the system with 300 particles (150 protons + 150 neutrons). The results 
obtained with the second method for the system with 512 particles (256 protons + 256 
neutrons) are totally consistent with those presented here and give the same equation of 
state with the same critical point. 

After the system has thermalized, we calculate the equation of state of our system by 
means of the virial theorem. The pressure of a system of classical particles interacting 
through two-body forces is given by [@,[7|,|§ 



where V is the volume of the box, p = A/V the density, Fij = —dVij/dr is the force 
acting between particles i and j, and n2(ri,r2) is the pair distribution function. E\ is the 
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a maxwellian velocity distribution by means of a Metropolis sampling ||. The system 





kinetic energy of particle i, and the brackets < • > denote the time average. In the above 
equation, the Boltzmann constant k is set equal to 1. To calculate the pressure, we let the 
system evolve after the thermalization for additional time and take the average over this 
time evolution. The density of the system is fixed, while the temperature is given in the 
equilibrium phase by 2/3 of the kinetic energy. 

In Fig. 3 we have plotted the pressure P versus the final temperature of the system T 
for several values of the density p ranging from 0.02 to 0.24 fm~ 3 . We see from the figure 
that for all the values of the density considered, the pressure P shows a linear behaviour 
versus temperature, 

P(p,T) = a(p)T + b(p) (3) 

Note that the temperatures considered here are that of the thermalized system (at equilib- 
rium). The values of the fitting parameters a(p) and b(p) are plotted in the upper and lower 
parts of Fig. 4, respectively, versus the density p (solid dots). At this point, to determine 
the equation of state, we fitted the parameters a(p) and b(p) with polynomials in p with the 
minimum number of fitting parameters and the following constraints: i) All the isotherms 
have a zero pressure at zero density P(p — 0, T) — 0, which means a(0) = 6(0) = 0; ii) 
The isotherm T = MeV has a zero slope at zero density f^-| T=0p=0 = 0, which means 
f^lp=o = 0- The resulting polynomial fits give: 

P(p, T) = (a lP + a 2 p 2 )T + b 2 p 2 + b 3 p 3 

= (0.96p + 7. 13p 2 )T - 87.0p 2 + 646.0p 3 (4) 

The fits to a(p) and b(p) are shown by solid lines in Fig. 4. Note that the parameter a x 
is almost equal to 1 which means that our equation of state is in agreement with that of a 
classical gas at high temperature and low density p — > 0,T — > oo =>- P ~ pT. Turning off 
the two-body potential, the law of ideal gases is respected within less than 1%. We have 
plotted in Fig. 5 the pressure P, eq.(4), versus density for several values of the temperature 
T. The dashed line in the figure delimits the isothermal spinodal region given by 
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The obtained equation of state eq.(|) has a critical point 

T c = 2.53 MeV ; p c = 0.04 fm^ ; P c = 0.029 MeV fm' 3 (6) 

Note that the ratio P c /T c p c = 0.32 is in good agreement with the experimental values for real 
gases which give for example 0.287 for CO2, 0.290 for Xe and 0.305 for 4 He. We have also 
plotted in Fig. 6 the isothermal spinodal line (eq.@) in the plane (p/p c , T/T c ) together 
with the experimental values for eight different substances in the gas-liquid coexistence 
region fLOf . The critical data for all these substances are quite varied (see Ref. 01), but the 
reduced data points fall on a universal curve. The right branch refers to the liquid phase, 
and the left branch to the gas phase, and both come together at the critical point. Our 
coexistence curve falls also, with some small differences, on the same universal curve. The 
dashed line shows Guggenheim's fit [|TTJ|] . 

As already stated, using the second method of initialization and thermalization for the 
system with 512 particles we obtained the same equation of state eq.@ with less that 5% 
difference with the first method using a system with 300 particles. This allows us to conclude 
that the equation of state of our system does not depend on the number of particles we have 
used in the box, and on the way how the system has reached equilibrium. 



III. MASS DISTRIBUTIONS 

In this section, our aim is to generate the mass distributions at the critical point and 
around it to see whether one gets a power law at the critical point as predicted by the 
droplet model of Fisher |3J], and as observed in CMD for an expanding classical system 
||. To do so, we generate many events near the critical temperature, and at the critical 
density p c = 0.04 fm~ 3 . Unfortunately our algorithm of cluster recognition which is of 
the percolation-type, is not able to identify the fragments inside the box, even at such 
small density [|ll|]. In this algorithm, the fragments are defined as follows. It is possible to 
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go from any particle i to another particle j belonging to the same fragment by successive 
interparticle jumps of a prescribed distance d or less (in these calculations, d is set equal 
to the range of the two-body potentials d = r c = 5.4 fm), but any path to a particle k 
belonging to another fragment contains one or more jumps of distance greater than d. This 
kind of cluster recognition algorithms is not indicated for identifying fragments confined in 
a box because the results are very dependent on the choice of the distance d. However, for 
expanding systems with fragments flying away from each other, these algorithms work well, 
particularly by following the fragments in time. To overcome this inherent problem, after 
the system has thermalized and the fragments are formed inside the box, we take off the 
box and let the system expand and the fragments fly away for enough time to allow a good 
identification of the fragments. 

Because this procedure is very CPU-time consuming, we have generated only 100 events 
of the system with 216 particles (108 protons + 108 neutrons) at each temperature and 
calculated the mass distributions from these events. The results of these calculations are 
drawn in Fig. 7, where we have plotted the mass distributions obtained from the expansion 
of system prepared at the critical density p c = 0.04 fm~ 3 and temperatures (at thermal- 
ization) T = 2.0 MeV, T c = 2.54 MeV and T = 3.92 MeV. At the critical temperature, 
apart for large masses where one observes large fluctuations due to the limited number of 
events generated, one observes a power law, consistent with the droplet model of Fisher, 
with a power r = —2.23, the same as the one found for the expanding finite system passing 
during its expansion through the critical point. Below the critical temperature, one ob- 
serves a mass distribution with a " U" shape characteristic to undercritical events where one 
sees some remaining of the initial system, and above the critical temperature one observes 
a rapidly decreasing mass distribution with an exponential shape characteristic to highly 
excited systems going to complete vaporization. 
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IV. FISHER'S DROPLET MODEL AND THE EQUATION OF STATE 



Fisher's droplet model has successfully been applied to describe the fragmentation of 
excited classical systems which expand and, depending on the excitation energy show dif- 
ferent dynamical evolutions, from evaporation like processes (small excitation energy) to 
multifragmentation and complete vaporization of the system (for large excitation energies) 

I- 

In the following, using the grand partition function Fisher has derived for a large system 
of individual particles in which clusters of two, three and more particles are bound together 
by the attractive forces, and in mutual statistical equilibrium with the remaining parts of 



the system and the surrounding gas [12]. In this model, the grand partition function 



z 



V,P) = Y^ z N Q N (V, 0) (7) 



N=0 



where z = eP^ is the fugacity, (3 = 1/T and Qn(V,P) is the canonical partition function, is 
given by 

JnS0j,V,T) = £ y (^) (8) 

A=l 

where Y(A) is the probability to find a cluster of size A in the system (cluster size distribu- 
tion) and is given by |12| 



Y(A) = Y A- T exp [fifiA - (3bA 2/3 ' 

= Y A- T B A2/3 S A (9) 

where // = \i g — fii is the difference between the chemical potentials of the gas and liquid 
phases, b denotes the surface contribution to the Gibbs free energy and is proportional to 
the surface tension. The term A~ T has been introduced by Fisher to take into account the 
fact that the droplet surface closes on itself. The parameter r is related to some critical 
exponents through scaling laws of critical phenomena and thus has a constant value in the 
range 2 < r < 2.5 P|,|i"2"|]. Knowing the grand partition function, the equation of state of the 
system is given by 
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ryi rp OO 

P = -ln~(»,V,T) = -J2Y(A) 

V V A=l 

< N > = z—lnE(v, V, T) = E AY ( A ) ( 10 ) 

° Z A=l 



and the isothermal compressibility is given by, 



VT d 2 



1 E A 2 Y{A) (11) 



T < N > 2 A=1 



Defining the moments of the cluster size distributions Y(A) by 



M k = E A k Y{A) (12) 

A=l 



gives for the previous quantities 



<N>=M 1 ^ P =^ = f (13) 



and 



XT = ^~ 2 M 2 (14) 
From eqs.(|13|), one obtains the equation of state in terms of the zeroth and first moments of 



the cluster size distribution 

P = TpM /M 1 (15) 

More generally, the moments of the cluster size distributions are related to Fisher's grand 
partition function by: 

M k = (r-^j lnS(ji,V,T) (16) 

The equation of state eq. fllBD is useless because of the complexity in estimating in a 
model independent way the cluster size distributions Y(A) (eq.@) and hence the moments 



Mfc. However, an estimate of the moments in the vicinity of the critical point is possible 
by making a few approximations. In the formula for the cluster size distributions eq.(^), 
for temperatures equal or larger than the critical temperature T > T c , the surface term b 
vanishes because it is proportional to the surface tension which is equal to for T >T C and 
the difference of chemical potentials \i is negative 



Y(A) = Y A~ T exp{(3pA) T>T C 



(17) 



We can now use the well known formula by Robinson derived for the Bose-Einstein integral 



HI : 



J2 e~ Aa A-° = T(l - a)a°- 1 + ]T 



A=l 



n=0 



III 



-((a-n)a r 



(18) 



where a > 0. 

Using this formula, the moments of the cluster size distribution become for T > T c : 



M k = Y 



T(l + k- r)\(3^ 1+k -^ + £ [ — 1 L C(t - k - n)\Pii\ n 



n=0 



111 



(19) 



The second term in the rhs of the above equation is regular and hence, while the zeroth and 
first moments are always finite because r is limited between 2 and 2.5, the higher moments 
diverge as 



M k oc \(3fi\ 



-(l+fc-T) 



(20) 



By assuming \i — (T — T c ) u , one recovers the formula Campi has derived for the moments 
at the critical point || 



M, oc IT _7T|-"(i+*-r) 



(21) 



In the vicinity of the critical point, the difference of chemical potentials [i tends to zero 
and one can consider the lowest orders only in the expansion of the moments M k in terms 
of powers of /x, eq.(|l9|), which gives for the pressure P and density p (assuming r = 2.33): 



P = q T [3.072|/3/i| 4/3 + 1.417 - 3.631|/fyt| + ... 



p = qo 



-4. 086 1 Pfj, | 1/3 + 3.631 + 0.966|/?//| + ... 



(22) 
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where go = Yq/V. At the critical point, p = and one gets, 



P c = lA17q T c 

p c = 3.631g (23) 



Subtracting eqs. (|22|) and (p23|) gives 



P-P c = q T [3.072|/^| 4/3 - 3.631\pp\ + ...] + lA17q (T - T t 



p- p c = qo 



-4.086|/3/i| 1/3 + 0.966|/3/x| + ...1 (24) 



Combining the above equations and taking the lowest orders in terms of powers of p, we get 
the following equation of state, valid in the vicinity of the critical point 



P-P c = 1.417g (T - T c ) + °-^T(p - p c f + °-^T(p - p c ) 4 (25) 



Using the numerical values of the critical temperature, density and pressure, obtained 
in section II for the system in a box with periodic boundary conditions (eq.@), we have 
plotted in Fig. 8, the critical isotherm of the above equation of state in the vicinity of the 
critical point (dashed line). In the same figure, we have also plotted the critical isotherm of 
the equation of state obtained in section II eq.(f|) (solid line). The critical point is indicated 
in the figure by the open circle. We see from the figure that the two curves go on each other 
in the vicinity of the critical point, while they start to deviate going away from the critical 
region. 

From eq . (|16D , we see that the isothermal compressibility xt is related to the second 
moment Mi- At the critical point, the isothermal compressibility becomes infinite, and so 
the second moment diverges, which means that the fluctuations of cluster sizes expressed by 
the second moment M 2 become very large at the critical point. Defining the relative variance 
72, introduced by Campi to have more insight in the shape of the cluster size distributions 
and to indicate the critical behaviour 0, 

M 2 M 

72 = ~w (26) 
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one expects that this quantity diverges at the critical point. Of course, if one considers a 
finite system, the moments remain finite, and instead of a divergence, the relative variance 
72 shows a peak around the critical point. To study finite size effects on the relative variance 
72, we calculate the moments of cluster size distributions Y(A) using for the parameters Yq, 
t, B and S appearing in eq. (|9[) the values obtained by fitting the fragment mass distributions 
obtained in Molecular Dynamics simulations of an expanding system with different initial 
excitation energies (temperatures) |§. In Fig. 9 we have plotted the relative variance 72 
versus the initial temperature (excitation energy) for different values of the total size of 
the system A to t- One can see the signature of the critical point for the initial temperature 
T = 5 MeV where one exactly obtains a power law in the fragment size distribution. One 
also sees that the reduced variance diverges more and more sharply as A tot increases. It is 
also interesting to compare the critical ratio P c /p c T c = Mq/Mi\t c calculated at the point 
where one observes the peak for 72, with other systems. For A to t = 100 and A to t = 10 8 , one 
gets Mq/Mi = 0.299 and 0.298, respectively. Our numerical simulation of the system in a 
box gives 0.317, and Van der Waals equation gives 0.375. Real gases give 0.287 for C0 2 , 
0.290 for Xe and 0.305 for 4 He [|HJ. 



V. CONCLUSIONS 

We have calculated by means of the virial theorem the equation of state of an infinite 
classical system in a cubic geometry with periodic boundary conditions and constant density. 
This equation of state has a critical temperature of 2.53 MeV and a critical density of 0.036 
fm~ 3 , about 1/3 the ground state density. We note that by using different numbers of 
particles in the box and different initialization and thermalization procedures, we obtain 
the same EOS with the same critical point. Another point is the decreasing of the critical 
temperature with the number of particles of the system. Several investigations of finite size 
effects on the critical temperature using Skyrme-type interactions |TJ[ have concluded that 



the critical temperature goes down when decreasing the number of particles of the system. In 
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our present work, it appears that this is not true and we have obtained for an infinite system 
(with periodic boundary conditions) the same critical temperature as the one observed in a 
previous work where we have studied, using the same two-body interaction the expansion 
of a finite system with 100 particles, passing during its evolution through the critical point. 
Furthermore, we have obtained at the critical point a fragment size distribution showing a 
power law with r = 2.23. The same power law was also observed for the expanding finite 
system for the critical evolution. 

Finally, by making use of the grand partition function of the droplet model of Fisher, we 
have studied the critical behaviour of the moments of fragment size distributions and found 
the same behaviour as the one calculated by Campi. We have also extracted the equation 
of state of this model around the critical point shown that it is in good agreement with that 
calculated for the infinite classical system. 
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FIGURES 

FIG. 1. Kinetic energy (upper part), potential energy (lower part, solid line) and total energy 
(lower part, dashed line) are plotted versus time for the system with 512 particles using the second 
method for thermalization (see text). The density is fixed at p = 0.16 fm~ 3 and the system is 
initialized at a temperature T = 1.66 MeV. 

FIG. 2. Same as Figure 1 but the density of the system is fixed at p = 0.04 fm~ 3 . 

FIG. 3. The pressure P of the system, calculated using the virial theorem, is plotted versus 
temperature T for several values of the density ranging from 0.03 to 0.24 fm~ 3 . The points indicate 
the results of the numerical simulations, and the solid lines the linear fits. The temperatures 
considered are the ones at equilibrium. 

FIG. 4. The values of the fitting parameters a(p) (upper part) and b(p) (lower part) are plotted 
versus density p (solid points). The solid lines show the polynomial fits to these values (see text). 

FIG. 5. The pressure P, calculated using the obtained equation of state eq.(4), is plotted versus 
density p for several values of temperature T. The dashed line indicates the isothermal spinodal 
line of the EOS. The critical isotherm is indicated by the small-dashes line. 

FIG. 6. Reduced temperature versus reduced density in the gas-liquid coexistence region for 
our system (solid line) and for eight different substances; Ne, A, Kr, Xe, N2, O2, CO and CH4 
(open circles). The dashed line shows Guggenheim's fit. 

FIG. 7. Mass distributions obtained for the fragmentation of the system with 216 particles with 
tree different temperatures (at equilibrium) T = 2.0 MeV (upper part), T c = 2.54 MeV (central 
part) and T = 3.92 MeV (lower part). 
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FIG. 8. Pressure versus density. The solid line shows the critical isotherm of the equation of 
state obtained in Section II by numerical simulations of a classical system, and the dashed line the 
critical isotherm of the analytical EOS obtained in Section IV using Fisher's droplet model. The 
open circle indicates the critical point. 

FIG. 9. The reduced variance 72 is plotted versus the initial temperature for the expansion of 
a finite system (see text) for several values of the size of the system Atot- 
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